Density functional method for nonequilibrium electron transport 
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We describe an ab initio metiiod for calculating the electronic structure, electronic transport, and 
forces acting on the atoms, for atomic scale systems connected to semi-infinite electrodes and with an 
applied voltage bias. Our method is based on the density functional theory (DFT) as implemented in 
the well tested Siesta approach (which uses non-local norm-conserving pseudopotentials to describe 
the effect of the core electrons, and linear combination of finite-range numerical atomic orbitals 
to describe the valence states). We fully deal with the atomistic structure of the whole system, 
treating both the contact and the electrodes on the same footing. The effect of the finite bias 
(including selfconsistency and the solution of the electrostatic problem) is taken into account using 
nonequilibrium Green's functions. We relate the nonequilibrium Green's function expressions to the 
more transparent scheme involving the scattering states. As an illustration, the method is applied 
to three systems where we are able to compare our results to earlier ab initio DFT calculations or 
experiments, and we point out differences between this method and existing schemes. The systems 
considered are: (1) single atom carbon wires connected to aluminum electrodes with extended or 
finite cross section, (2) single atom gold wires, and finally (3) large carbon nanotube systems with 
point defects. 
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I. INTRODUCTION 

Electronic structure calculations are today an impor- 
tant tool for investigating tho, physics and chemistry of 
new molecules and materials.ld An important factor for 
the success of these techniques is the development of first 
principles methods that make reliable modeling of a wide 
range of systems possible without introducing system de- 
pendent parameters. Most methods are, however, limited 
in two aspects: (1) the geometry is restricted to either 
finite or periodic systems, and (2) the electronic system 
must be in equilibrium. In order to address theoreti- 
cally the situation where an atomic/molecular-scale sys- 
tem (contact) is connected to bulk electrodes of requires 
a method capable of treating an infinite and non-periodic 
system. In the case where a finite voltage bias applied to 
the electrodes drives a current through the contact, the 
electronic subsystem is not in thermal equilibrium and 
the model must be able to describe this nonequilibrium 
situation. The aim of the present work is to develop a 
new first principles nonequilibrium electronic structure 
method for modeling a nanostructure coupled to exter- 
nal electrodes with diff'erent electrochemical potentials 
(we will interchange the terms electrochemical potential 
and Fermi fe?;e/ throughout the paper). Besides, we wish 
to treat the whole system (contact and electrodes) on the 
same footing, describing the electronic structure of both 
at the same level. 

Our method is based on the Density Functional Theory 
(DFT) .ElElu'B In principle the exact electronic density and 
total energy can be obtained within the DFT if the exact 
exchange-correlation (XC) functional was available. This 
is not the case and the XC functional has to be substi- 
tuted by an approximate functional. The most simple 
form is the Local Density Approximation (LDA), but re- 
cently a number of other more complicated functionals 



have been proposed, which have been shown to genei^ 
ally improve the description of systems in equilibrium. Q 
There is no rigorous theory of the validity range of these 
functionals and in practice it is determined by testing the 
functional for a wide range of systems where the theoret- 
ical results can be compared with reliable experimental 
data or with other more precise calculations. 

Here we will take this pragmatic approach one step 
further: We will use not only the total electron density, 
but the Kohn-Sham wave functions as bona fide single- 
particle wave functions when calculating the electronic 
current. Thus we assume that the commonly used XC 
functionals are able to describe the electrons in non- 
equilibrium situations where a cutrent flow is present, 
as in the systems we wish to study.H This mean-field-like, 
one-electron approach is not able to describe pronounced 
many-body effects which may appear in in some cases 
during the transport process. Inelastic scattering e.g. by 
phononaa will not be considered, either. 

Except for the approximations inherent in the DFT, 
the XC functional, and the use of the Kohn-Sham wave 
functions to obtain a current, all other approximations in 
the method are controllable, in the sense that they can 
be systematically improved to check for convergence to- 
wards the exact result (within the given XC functional). 
Examples of this are the size and extent of the basis set 
(which can be increased to completeness), the numeri- 
cal integration cutoffs (which can be improved to conver- 
gence) , or the size of the electrode buffer regions included 
in the selfconsistent calculation (see below). This mean- 
field-like, one-electron approach is not able to describe 
pronounced many-body effects which may appear in in 
some cases during thcrtransport process. Inelastic scat- 
tering e.g. by phononaa will not be considered, either. 

Previous calculations for open systems have 
in most cases been based on semi-empirical 
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approaches .100111110000300011 The first 
non-equihbrium calculations with a full self-consistent 
DFT description of the entire system have— eaaployed 
the jellium approximation in the electrodes Other 
approaches have used an equilibrium first principles 
Hamiltonian for the nanostructure and described the 
electrodes by inrlnHirig soTTji-pmpiriral sclfcnergies on 
the outermost atoms.L§LjO Lately, there have been 
several approaches which treat the .''^^Tfj system on 
the same footing, at the atomic levelcjO'Ej but so far 
only one of the approaches have been applied to the 
nonequilibrium situation where the-sxternal leads have 
different electrochemical potentials. lJO 

The starting point for our implcpaentation is the 
Siesta electronic structure approach^ In this method 
the effect of the core electroas,is described by soft norm- 
conserving pseudopotentialsEJ and the electronic struc- 
ture of the valence electrons is expanded in a, bjfisis set 
of numerical atomic orbitals with finite rangeOcJ The 
quality of the basis set can be improved at will Jay using 
multiple-C orbitals, polarization functions, etc.,E3 allow- 
ing to achieve convergence of the results to the desired 
level of accuracy. Siesta has been tested, ta a wide va- 
riety of systems, with excellent resultsO'cJ The great 
advantage of using (pbitals with finite range (besides the 
numerical efficienc};^) , is that the Hamiltonian interac- 
tions are strictly zero beyond some distance, which al- 
lows to partition the system unambiguously, and define 
regions where we will do different parts of the calculation 
as we describe in the Sections II-|rv|. Besides, the Hamil- 
tonian takes the same form as in empirical tight-binding 
calculations, and therefore the techniques developed in 
this context can be straightforwardly applied. 

We have extended the Siesta computational pack- 
age to nonequilibrium systems by calculating the den- 
sity matr ix ,, with la nonequilibrium Greens functions 
technique J130'l3e3 We have named this nonequilibrium 
electronic structure code TranSiesta. Preliminary 
results obtained with TranSiesta were presented in 
Ref. 41, Here we give a detailed account of the techni- 



cal implementation and present results for the transport 
properties of different atomic scale systems. One of the 
authors (JT) has been involved iur-tiie independent de- 
velopment of a package, McDCALJ23 which is based on 
similar principles, but with some differences in implemen- 
tation. We compare results obtained with the two pack- 
ages for a carbon wire connected to aluminum electrodes 
and show that they yield similar results. We present new 
results for atomic gold wire systems which are one of 
the most studied atomic scale conductors, and finally we 
present results for transport in nanotubes with defects. 

The organization of the paper is the following. In the 
first part of the paper we describe how we divide our 
system into the contact and electrode parts and how we 
obtain the density matrix for the nonequilibrium situa- 
tion using Green's functions. Here we also discuss the 
relation between the scattering state approach and the 
nonequilibrium Green's function expression for the den- 



sity matrix. Then we describe how this is implemented 
in the numerical procedures and how we solve the Pois- 
son equation in the case of finite bias. In the second 
part of the paper we turn to the applications where our 
aim is to illustrate the method and show some of its ca- 
pabilities rather than presenting detailed analysis of our 
findings. We compare our results with other ab initio cal- 
culations or experiments for (1) carbon wires connected 
to aluminum electrodes, (2) gold wires connected to gold 
electrodes, and finally (3) infinite carbon nanotubes con- 
taining defects. 



II. SYSTEM SETUP 

We will consider the situation sketched in Fig. l|a. Two 
semi-infinite electrodes, left and right, are coupled via a 
contact region. All matrix elements of the Hamiltonian 
or overlap integrals between orbitals on atoms situated 
in different electrodes are zero so the coupling between 
the left and right electrodes takes place via the contact 
region only. 




(a) 



Semi-infinite bulk 



(b) 




FIG. 1: (a) We model the Contact (C) region coupled to two 
semi-infinite Left (L) and Right {R) electrodes. The direction 
of transport is denoted by z. (b) We only describe a finite 
section of the infinite system: Inside the L and R parts the 
Hamiltonian matrix elements have bulk electrode values. The 
external {Buffer) region, B, is not directly relevant for the 
calculation. 

The region of interest is thus separated into three parts. 
Left (L), Contact (C) and Right (R). The atoms in L 
(R) are assumed to be the parts of the left (right) semi- 
infinite bulk electrodes with which the atoms in region 
C interact. The Hamiltonian is assumed to be converged 
to the bulk values in region L and R along with the den- 
sity matrix. Thus the Hamiltonian, density and overlap 
matrices only differ from bulk values in the C, C — L and 
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C — R parts. We can test this assumption by including 
a larger fraction of the electrodes in C (so the L and R 
regions are positioned further away from the surfaces in 
Fig.0). 

In order to obtain the transport properties of the sys- 
tem, we only need to describe the finite L — C—R part of 
the infinite system as illustrated in Fig. |l|b. The density 
matrix which describes the distribution of electrons can 
be obtained from a series of Green's function matrices 
of the infinite system as we will discuss in detail in Sec- 
tion III . In principle the Green's function matrix involves 
the inversion of an infinite matrix corresponding to the 
infinite system with all parts of the electrodes included. 
We are, however, only interested in the finite L — C — R 
part of the density matrix and thus of the Green's func- 
tion matrix. We can obtain this part by inverting the 
finite matrix. 
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where H/,, H^j and Hp are the Hamiltonian matrices in 
the L, R and C regions, respectively, and (Vr) is 
the interaction between the L (R) and C regions. The 
coupling of L and R to the remaining part of the semi- 
infinite electrodes is fully taken into account by the self- 
energies. El and Tir. We note that to determine V^, 
Vij and He, we do not need to know the correct density 
matrix outside the L — C — R region, as long as this does 
not influence the electrostatic potential inside the region. 
This is the case for metallic electrodes, if the L — C — R 
region is defined sufficiently large so that all the screening 
takes place inside of it. 

The upper and lower part of the Hamiltonian (HL(i?,) + 
^L(ii)) are determined from two separate calculations for 
the bulk systems corresponding to the bulk of the left 
and right electrode systems. These systems have periodic 
boundary conditions in the z directions, and are solved 
using Bloch's theorem. From these calculations we also 
determine the self-energies by cutting the electrode sys- 
tems into two, semi-infinite pieces using either the ideal 
constructioncS or the efficient recursion method.c3 

The remaining parts of the Hamiltonian, V^, Vj^ and 
H(7, depend on the non-equilibrium electron density and 
are determined through a selfconsistent procedure. In 
Section III we will describe how the non-equilibrium den- 



sity matrix can be calculated given these parts of the 
Hamiltonian, while in Section IV we show how the effec- 



tive potential and thereby the Hamiltonian matrix ele- 
ments are calculated from the density matrix. 



III. NON-EQUILIBRIUM DENSITY MATRIX 

In this section we will first present the relationship 
between the scattering state approach and the non- 
equilibrium Green's function expression for the non- 



equilibrium electron density corresponding to the situ- 
ation when the electrodes have different electrochemical 
potentials. The scattering state approach is quite trans- 
parent and has been used for non-equilibipim first princi- 
ples calc.ulatinjis by McCann and Brownjlj Lang and co- 
workers,L3oc3 and Tsukada and co-workers Oc3a All 
these calculations have been for the case of model jel- 
lium electrodes and it is not straightforward how to ex- 
tend these methods to the case of electrodes with a real- 
istic atomic structure and a more complicated electronic 
structure or when localized states are present inside the 
contact region. The use of the non-equilibrium Green's 
functions combined with a localized basis set is able to 
deal with these points more easily. 

Here we will start with the scattering state approach 
and make the connection to the non-equilibrium Green's 
function expressions for the density matrix. Consider the 
scattering states starting in the left electrode. These are 
generated from the unperturbed incoming states (labeled 
by I) of the uncoupled, semi-infinite electrode, ■0/', using 
the retarded Green's function, G, of the coupled system, 

M^) = i^ii^) + fdy G(x, y) VUy) 0,°(y) . (2) 



As in the previous section there is no direct interaction 
between the electrodes: 



V{r) 
Vi?(r)V^?(rl 



^L(r)^°(f) = 0. 



(3) 
(4) 



Our non-equilibrium situation is described by the fol- 
lowing scenario: The states starting deep in the left /right 
electrode are filled up to the electrochemical potential of 
the left (right) electrode, /il (m_r). We construct the den- 
sity matrix from the (incoming) scattering states of the 
left and right electrode: 

D{x,y) = ^ipi{x)tlJi{y)nF{ei - hl) 
I 

+ ^ipr-{x)ip*{y)nF{er- - fiR) , (5) 

r 

where index / and r run over all scattering states in the 
left and right electrode, respectively. Note that this den- 
sity matrix only describes states in C which couple to 
the continuum of electrode states - we shall later in Sec- 



tion III D return to the states localized in C. 



A. Localized non-orthogonal basis 

Here we will rather consider the density matrix de- 
fined in terms of coefficients of the scattering states 
with respect to the given basis (denoted below by Greek 
subindexes) 



(6) 
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Thus (|) and (|) read, 

= < + E (G(^)V)^, , z = ei + i5, (7) 



(8) 



The basis is in general non-orthogonal but this will not 
introduce any further complications. As for the Hamilto- 
nian, we assume that the matrix elements of the overlap, 
S, are zero between basis functions in L and R. The 
overlap is handled by defining the Green's function ma- 
trix G(z) as the inverse of (zS — H), and including the 
term — zS in the perturbation matrix V. To see this we 
use the following equations, 



[ejSo-Ho]cO = 0, 
hS-H]Q = 0, 
[zS-H]G(2) = 1. 

With these definitions we see that (0) is fulfilled, 

[eiS-U] ci - [e,S-H] cO + Vc? = 0, 

when 



V = H Hn 



ei(S-So). 



(9) 
(10) 
(11) 

(12) 
(13) 



The use of a non-orthogonal basis is described further in 



Refs. m and 49 



The density matrix naturally splits into left and right 
parts. The derivations for left and right are similar, so we 
will concentrate on left only. It is convenient to introduce 
the left spectral density matrix, p^. 



^cif,c*i^5{e - ei), 



(14) 



and likewise a right spectral matrix pR. The density 
matrix is then written as, 

/oo 
d£p^^{e)nF{e - Ml) + p^'Je)nF{e - ^ir) . 
-OO 

(15) 

As always we wish to express in terms of known 
(unperturbed) quantities, i.e., c^^, and for this we use 
Equation (^. Since we are only interested in the den- 
sity matrix part corresponding to the scattering region 
[L — C — R), we note that the coefficients for the un- 
perturbed states are zero for basis functions (/i) within 
this region. Thus 



(16) 



where v is inside the bulk of the left electrode. Inserting 
this in Eq. (|l|) we get. 



p^,(£)= G(e)-Im[Vg^(e)Vt] G\e, 



(17) 



Here we use the unperturbed left retarded Green's func- 
tion. 



„0 „0* 

e — ei + i6 



and the relation 

[g^(e) - (g^(e))t]^^ ^ 2mY,clct 5{e - si) . 

and that g = g'^ due to time-reversal symmetry. 
We can identify the retarded self-energy, 

Ei(e) ^ [Vg^(e)Vt] , 

TUz) = z[SL(e)-EL(e)t]/2, 

and finally we express as, 

P^.(e) = -(G(e)rL(£)Gt(£)) 



(18) 



(19) 



(20) 
(21) 



(22) 



and a similar expression for p^. Note that the S, F and G 
matrices in the equations above are all matrices defined 
only in the scattering region L — C ~ R which is desirable 
from a practical point of view. The G matrix is obtained 
by inverting the matrix in Eq. 0. 

The expression derived from the scattering states is the 
same as one would get from a non-equilibrium Green's 
function derivation, see e.g. Refs. |3^, where D is ex- 
pressed via the "lesser" Green's function. 
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(23) 



which includes the information about the non- 
equilibrium occupation. 



B. Complex contour for the equilibrium density 
matrix 

In equilibrium we can combine the left and right parts 
in Eq. (H), 



GEG^ = -G [E-Stl Gt 
2 ^ 



-G [(G)-i - (Gt)-i] Gt 



= -Im[G] 



(24) 



where S includes both E/, and S^, and time-reversal 
symmetry (G^ = G*) was invoked. With this Eq. (|l|) 
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reduces to the well-known expression, 



D 



1 



TT 



de Im[G(£ + iS)] npi^s — A^) 

deG{e + i5)nF{e- ^Ji)\. (25) 



Im 

TT 



The invocation of time-reversal symmetry makes D a real 
matrix since D* = D-^ = D. 

At this point it is important to note that we have ne- 
glected the infinitesimal i6 in Eq. (p4[). This means that 
the equality in Eq. (|2j) is actually not true when there 
are states present in C which do not couple to any of 
the electrodes, and thus = Tn = and pl — PR — 
for elements involving strictly localized states. The lo- 
calized states cannot be reached starting from scattering 
states and are therefore not included in Eq. (^5|), while 
they are present in Eq. (25). We return to this point in 
Section HID below. 




FIG. 2: The closed contour: L (]oo + iA;EF - j + iA[), C, 
and [EB + iS;<x + iS] enclosing the Fermi poles (black dots). 

All poles of the retarded Green's function G(z) are ly- 
ing on the real axis and the function is analytic otherwise. 
Instead of doing the integral in Eq. ( p5| ) (corresponding 
to the dotted line in Fig. ||), we consider the contour 
in the complex plane defined for a given finite tempera- 
ture shown by the solid line in Fig. g. Indeed, the closed 
contour beginning with line segment L, followed by the 
circle segment C, and running along the real axis from 
{EB + iS) to (oo -I- iS), where EB is below the bottom 
valence band edge, will only enclose the poles of np{z) 
located at Zi, = i{2v + \)nkT. According to the residue 
theorem, 



^ dz G{z) np{z ^ p) = — 27rj kT ^ ^ G{ziy) . 



(26) 



where we use that the residues oi np are —kT. Thus, 

J^gdeG{e + i5)np{e- p) (27) 
= — Jq^i^ dz G{z) np{z — p) — 2ni kTJ2z„ G(z^) . 

The contour integral can be computed numerically for 
a given finite temperature by choosing the number of 
Fermi poles to enclose. This insures that the complex 
contour stays away from the real axis (the part close to 
EB is not important). The Green's function will behave 



smoothly sufficiently away from the real axis, and we 
can do the contour integral by Gaussian quadrature with 
just a minimum number of points, see Fig. |. The mam 
variation on L comes from np and it is advantageous to 
use ni^- as a weightfunction in the Gaussian quadrature.E2l 




FIG. 3: Typical points for Gaussian quadrature on the con- 
tour. On L we employ a quadrature with a weight function 
equal to the Fermi function. 



Numerical procedure for obtaining the 
non-equilibrium density matrix 



In non-equilibrium the density matrix is given by 



D 
D 



L 



(28) 



= Im 

TT 

• oo 



deG{e + i6)np{e - pl) , (29) 



EB 
R 



- , dep^^{e)[np{e- pR)-np{e- pL)]{;iO) 



or equivalently 
D - 



= Im 



EB 



(31) 

deG{e + i5)np{e — pn) , (32) 



de Pnu{^)[np{e - pl) - npis - Mi?)](33) 



The spectral density matrices, p^ and p^ , are not an- 
alytical. Thus only the "equilibrium" part of the den- 
sity matrix, D-^(D^), can be obtained using the complex 
contour. Furthermore, this is a real quantity due to the 
time-reversal symmetry, whereas the "non-equilibrium" 
part, A^(A^) cannot be made real since the scattering 
states by construction break time-reversal symmetry due 
to their boundary conditions. The imaginary part of Ai 
(A^) is in fact directly related to the local current.EiJ 
However, if we are interested only in the electron density 
and if we employ a basis-set with real basis functions (0^) 
we can neglect the imaginary part of D , 



D, 



(34) 
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To obtain (A^) the integral must be evaluated for 
a finite level broadening, iS, and on a fine grid. Even 
for small voltages we find that this integral can be prob- 
lematic, and care must be taken to ensure convergence 
in the level broadening and number of grid points. Since 
we have two similar expressions for the density matrix 
we can get the integration error from 

e^. - D« + A^, - (D;:, + A« ). (35) 

The integration error arises mainly from the real axis 
integrals, and depending on which entry of the density 
matrix we are considering either A^ or A^ can domi- 
nate the error. Thus with respect to the numerical im- 
plementation the two formulas Eq. (|2^) and (|3^) are not 
equivalent. We will calculate the density matrix as a 
weighted sum of the two integrals 

D^, = ^,,.(D;;,-|-A^^J + (1-u;^.)(D;?, + A^j356) 

(A^ )2 
- (A^J2 + (A«)2- 

The choice of weights can be rationalized by the following 
argument. Assume that the result of the numerical inte- 
gration is given by a stochastic variable A^ with mean 
value A^ and the standard deviation is proportional to 
the overall size of the integral, i.e. Var{A^) cx (A^)^. 
A numerical calculation with weighted integrals as in 
Eq. ( p6| ) will then be a stochastic variable with the vari- 
ance 

Var{i))(xw^{A^f + {l-wf{A^f. (38) 

The value of w which minimize the variance is the weight 
factor we use in Eq. (p7|). 



D. Localized states 

As mentioned earlier the signature of a localized state 
at eo in the scattering region is that the matrix elements 
of TL{eo),^R{£o) are zero for that particular state. Lo- 
calized states most commonly arise when the atoms in 
C have energy levels below the bandwidth of the leads. 
The localized states give rise to a pole at eo in the Green's 
function. As long as Eq < {^L,fJ'R} the pole will be en- 
closed in the complex contours and therefore included in 
the occupied states. If on the other hand the bound state 
has an energy within the bias window, i.e. fi^ < Bq < j-iji 
the bound state will not be included in the real axis in- 
tegral (A^, A-'^) and in the complex contour for D^, but 
it will be included in the complex contour for D^. Such 
a bound state will only be correctly described by the 
present formalism if additional information on its filling 
is supplied. These situations are rare and seldom encoun- 
tered in practice. 



IV. THE NONEQUILIBRIUM EFFECTIVE 
POTENTIAL 

The DFT effective potential consists of three parts: 
a pseudo potential T^s, the exchange correlation poten- 
tial Vxc, and the Hartree potential Vh- For Vps we use 
norm conserving TrouUier-Martins pseudopotentials, de- 
termined from standard procedures.Ej For Vxc we use the 
LDA as parametrized in Ref. ^ 

A. The Hartree potential 

The Hartree potential is a non-local function of the 
electron density, and it is determined through the Pois- 
son's equation (in Hartree atomic units), 

VVH(r) = -47rn(r) . (39) 

Specifying the electron density only in the C region of 
Fig. |l| makes the Hartree potential of this region unde- 
termined up to a linear ternS, 

VH{r)^(j){r) + d-r + h, (40) 

where 0(r) is a solution to Poisson's equation in region C 
and a , h are parameters that must be determined from 
the boundary conditions to the Poisson's equation. In the 
directions perpendicular to the transport direction {x,y) 
we will use periodic boundary conditions which fix the 
values of ax, ay. The remaining two parameters , b are 
determined by the value of the electrostatic potential at 
the L — C and C—R boundaries. The electrostatic poten- 
tial in the L and R regions could be determined from the 
separate bulk calculations, and shifted relative to each 
other by the bias V . With these boundary conditions 
the Hartree potential in the contact is uniquely defped, 
and could be computed using a real-space techniqueO or 
an iterative method.O 

However, in the present work, we have solved the Pois- 
son's equation using a Fast Fourier Transform (EFT) 
technique. We set up a supercell with the L — C — i? re- 
gion, which can contain some extra layers of buffer bulk 
atoms and, possibly, vacuum (specially if the two elec- 
trodes are not of the same nature, otherwise the L and 
R are periodically matched in the z direction). We note 
in passing that this is done so that the potential at the 
L — C and C—R boundaries reproduces the bulk values, 
crucial for our method to be consistent. For a given bias, 
V , the L and R electrode electrostatic potentials need to 
be shifted by V/2 and — F/2, respectively, and Vh will 
therefore have a discontinuity at the cell boundary. The 
electrostatic potential of the supercell is now decomposed 
as: 

VHir) = 0(f) - - 0.5). (41) 

where ^(r) is a periodic solution of the Poisson's equa- 
tion injtiie supercell, and therefore can be obtained using 
FFT's.eI 
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10 20 30 40 



z(a„) 

FIG. 4: (a) The induced external potential for slab calculation 
(full line), and in the TranSiesta calculation (dashed line). 
In the slab calculation the jump in external potential is in the 
middle of the vacuum region. The total potential (arbitrary 
units) is shown for reference (dotted line), (b) Induced den- 
sity. Potential and density is averaged in the surface plane. 
The density corresponds to one surface unit cell. 



sity and the electrostatic potential to be correct at the 
L ~ C and C — R boundaries, the and parts of 
the Hamiltonian (see Eq. (^) will not be correct. We 
therefore substitute and H/j with the Hamiltonian 
obtained from the calculation of the separate bulk elec- 
trode systems. Here it is important to note that the 
effective potential within the bulk electrode calculations 
usually are shifted rigidly relative to the effective poten- 
tial in the L and R regions, due to the choice of the 
parameter b in Eq. (^). However, the bulk electrode 
Hamiltonians and Hr can easily be shifted, using 
the fact that the electrode Fermi level should be similar 
to the Fermi level of the initial Siesta calculation for the 
BLCRB supercell. 

The discontinuity of the Hartree potential at the cell 
boundary has no consequence in the calculation: the 
Hamiltonian matrix elements inside the L — C — R region 
are unaffected because of the finite range of the atomic 
orbitals, and the Hamiltonian matrix elements outside 
the L — C — R region which do feel the discontinuity are 
replaced by bulk values (shifted according to the bias) . 

V. CONDUCTANCE FORMULAS 



To test the method we have calculated the induced 
density and potential on a "capacitor" consisting of two 
gold (111) surfaces separated by a 12 bohr wide tunnel- 
gap and with a voltage drop of 2 Volt. We have calcu- 
lated the charge density and the potential in this system 
in two different ways. First, we apply the present for- 
mulation (implemented in TranSiesta), where the sys- 
tem consists of two semi-infinite gold electrodes, and the 
Hartree potential is computed as described above. Then, 
we calculate a similar system, but with a slab geometry, 
computing the Hartree potential with Siesta, adding the 
external potential as a ramp with a discontinuity in the 
vacuum region. Fig. ^ shows the comparison of the re- 
sults for the average induced density and potential along 
the z axis. Since the tunnel gap is so wide that there 
is no current running, the two methods should give very 
similar results, as we indeed can observe in the Figure. 
We can also observe that the potential ramp is very effec- 
tively screened inside the material, so that the potential 
is essentially equal to the bulk one, except for the surface 
layer. This justifies our approach for the partition of the 
system, the solution of the Poisson's equation, and the 
use of the bulk Hamiltonian matrix elements fo the L and 
R regions (see below). 



B. The Hamiltonian matrix elements 

Having determined the effective potential we calculate 
the Hamiltonian matrix elements as in standard Siesta 
calculations. However, since we only require the den- 



Using the non-equilibrium Green's function formalism 
(see e.g. Refs. p^[|39|j40| and references therein) the cur- 
rent, /, through the contact can be derived 

/oo 
de {np{e - ^l) - npie - ^J-b.)) 
-OO 

Tr[TUe)G\e)Tn{e)G{e)] , (42) 

where Go = 2e^ /h. We note that this expression is-iuot 
general but is valid for mean-field theory like DFT.ES 
equivalent formula has been derived by Todorov et alrj 
(the equivalence can be derived using Eq. (|2l|) and the 
cyclic invariance of the trace in Eq. ^2|) . 

With the identification_pf the (left-to-right) transmis- 
sion amplitude matrix t,E3 

tie) = {Tnie)y/'G{e){rL{e)y/\ (43) 

( ^ ) is |Seen to be equivalent to the Landauer-Biittiker 
formulaEj for the conductance, G = I/V, 

G 1"°^ 

G{V) = — ^ / de[nF{e - hl) - npie - fin)] 

^ J —oo 

Tr [ttt] (e) , (44) 

The eigenchannels are defined -ib terms of the (left-to- 
right) transmission matrix tJ^iJ'Ell 

t = U^diag{|T„|}U[, (45) 

and split the total transmission into individual channel 
contributions, 

TTot = El^"l'- (46) 
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The collection of the individual channel transmissions 
{|t„P} gives a more detailed description of the con- 



duct an 
results 



Ml 



is useful for the interpretation of the 



VI. APPLICATIONS 



A. Carbon wires/ Aluminum (100) electrodes 

Short monoatomic carbon wires coupled to metal- 
lic electMidcs have recently been studied by Lang and 
AvourisE3E2l and Larade et ali^. Lang and Avouris used 
the Jellium approximation for the electrodes, while La- 
rade et al. used Al electrodes with a finite cross section 
oriented along the (100) direction. In this section, we 
will compare the TranSiesta method with these other 
first principles electron transport methods by studying 
the transmission through a 7-atom carbon chain coupled 
to Al(lOO) electrodes with finite cross sections as well as 
to the fuU Al(lOO) surface. 

We consider two systems, denoted A and B, shown in 
Fig. ^,b. System A consists of a 7-atom carbon chain 
coupled to two electrodes of finite cross section oriented 
along the Al(lOO) direction (see Fig. ||a). The electrode 
unit cell consists of 9 Al atoms repeated to z = ±oo. The 
ends of the carbon chain are positioned in the Al(lOO) 
hollow site and the distance between the ends of the car- 
bon chain and the first plane of Al atoms is fixed to be 
d — 1.0 A. In system B the carbon chain is coupled to 
two Al(100)-(2\/2 X 2\/2) surfaces with an Al-C couphng 
similar to system A. In this case the electrode unit cell 
contains 2 layers each with 8 atoms. For both systems 
the contact region (C) includes 3 layers of atoms in the 
left electrode and 4 layers of the right electrode. We 
use single-C basis sets for both C and Al to be able to 
compare with the resuLts, from McDCAL, which were ob- 
tained with that basis .Eil 

The conductance of system A is dominated by the 
alignment of the LUMO state of the isolated chain to the. 
Fermi level of the electrodes through charge transfer.Efl 
The coupling of the LUMO, charge transfer, and to- 
tal conductance can be varied ccptinuously by adjust- 
ing the electrode-chain separation.!^ For our value of the 
electrode-chain separation we get a charge transfer of 1.43 
and 1.28 e to the carbon wire in systems A and B, respec- 
tively. This is slightly larger than the values obtained by 
Lang and Avourialj for Jellium electrcuies, but in good 
agreement with results from McDCAL.Efl 

To facilitate a more direct comparison between the 
methods we show in Fig. ^ the transmission coefficient 
of system A calculated both within TranSiesta (solid) 
and McDCAL (dotted). For both methods, we have used 
identical basis sets and pseudopotentials. However, sev- 
eral technical details in the implementations differ and 
may lead to small differences in the transmission spec- 
tra. The main implementation differences between the 
two methods are related to the calculation of Hamilto- 



nian parameters for the electrode region, the solution of 
the Poisson's equation, and-t|he complex contours used to 
obtain the electron charge.E3 Thus, there are many tech- 
nical differences in the two methods, and we therefore 
find the close agreement in Fig. ^ very satisfactory. 

In Fig. [t] we show the corresponding transmission co- 
efficients for system B. It can be seen that the trans- 
mission coefficient for zero bias at e = /i is close to 1 
for both systems, thus they have similar conductance. 
However, the details in the transmission spectra differs 
much from system A. In order to get some insight into 
the origin of the different features we have projected the 
selfconsistent Hamiltonian onto the carbon orbitals, and 
diagonalised this subspace Hamiltonian to find the posi- 
tion of the carbon eigenstates in the presence of the Al 
electrodes. Within the energy window shown in Fig. ^ 
and ^ we find 4 doubly degenerate 7r-states (Stt, in, 5n, 
Gtt). The positions of the eigenstates are indicated above 
the transmission curves. Each doubly degenerate state 
can contribute to the transmission with 2 at most. Gen- 
erally, the position of the carbon tt states give rise to a 
slow variation in the transmission coefficient, and the fast 
variation is related to the coupling between different scat- 
tering states in the electrodes and the carbon 7r-states. 
For instance in system A, there are two energy intervals [- 
1.9,-1.7] and [0.7,1.4], where the transmission coefficient 
is zero, and the scattering states in these energy inter- 
vals are therefore not coupling to the carbon wire. Note 
how these zero transmission intervals are doubled at fi- 
nite bias, since the scattering states of the left and right 
electrode are now displaced. 

The energy dependence of the transmission coefficient 
is quite different in system B compared to system A. This 
is mainly due to the differences in the electronic struc- 
ture of the surface compared to the finite sized electrode. 
However, the electronic states of the carbon wires are 
also slightly different. We find that the 7r-states lie 0.2 
eV higher in energy in system B relative to system A. As 
mentioned previously, the charge transfer to the carbon 
chain is different in the two systems. The origin of this is 
related to a larger work function (~ leV) of the surface 
relative to the lead. We note that the calculated work 
function of the surface is 0.76 eV higher than the experi- 
mental workfunction of Al (4.4 eV),E3 which may be due 
to the use of a single-^ basis sst-jand the approximate 
exchange-correlation description.EH 

In Fig. |d, we show the changes in the effective poten- 
tial when a 1 V bias is applied. We find that the potential 
does not drop continuously across the wire. In system A, 
the main potential drop is at the interface between the 
carbon wire and the right electrode, while in system B the 
potential drop takes place at the interface to the left elec- 
trode. This should be compared to the Jellium results, 
where there_is a more continuous voltage drop through 
the system.Cj We do not yet understand the details of the 
origin of these voltage drops. However, it seems that the 
voltage drop is very sensitive to the electronic structure 
of the electrodes. Thus, we find it is qualitatively and 
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FIG. 6: (a) Zero bias transmission coeSicient, r(i5,0V'), for 
the 7-atom carbon chain with finite cross section Al(lOO) 
electrodes (system A), (b) Transmission coefficient at 1 V, 
T{E,IV). Solid lines show results obtained with TranSi- 
ESTA, and dotted lines results obtained with McDCAL. The 
vertical dashed lines indicated the window between hl and 
fiR. The position of the eigenstates of the carbon wire sub- 
system are also indicated at the top axis. 



FIG. 5: (a) The 7-atom carbon chain with finite cross sec- 
tion Al(lOO) electrodes (system A), (b) The carbon chain with 
Al(100)-(2\/2 X 2^2) electrodes(system B). (c) The effective 
potential of system A(dashed) and system B(solid), together 
with the effective potential of the corresponding bare elec- 
trode systems, (d) The selfconsistent effective potential for 
an external bias of 1 V (the zero bias effective potential has 
been subtracted). 



quantitatively important to have a good description of 
the electronic structure of the electrodes. 




-2-10123 
E-iii^+ii^)/2 (eV) 



B. Gold wires/Gold (111) electrodes 

The conductance of single atom gold wires is a 
benchmark in. atomic scale conduction. Since the first 
experimenta^aSEj numerous detailed studies of their 
conductance have been carried out through the 1990s 
until now (see e.g. Ref. for j^tjwadew). More re- 
cently the non-linear conductanceOJSEjyij has been in- 
vestigated and the atomic structureLallZlta of these sys- 
tems has been elucidated. Experiments show that chains 
containing more than 5 gold atomglI3 can be pulled and 
that these can remain stable for an extended period of 
time at low temperature. A large number of experiments 
employing different techniques and under a variety of 
conditions (ambient pressure and UHV, room and liq- 
uid He temperature) all show that the conductance at 
low bias is very close to 1 Go and several experiments 
point to the facLiJiat this is due to a single conductance 
eigenchannel P^P^tr^ 

Several theoretical inv c)s| . ifflat i| 0rasj 1pave addressed the 
stability and morphologyHEaEaoOEa and the conduc- 
tanceE3E3'E3 of atomic gold chains and contacts using 



FIG. 7: Transmission coefficients, T{E, OV) and T{E, IV) 
for the 7-atom carbon chain with A1(100)-(2a/2 x 2\/2) elec- 
trodes(system B). The position of the eigenstates of the car- 
bon wire subsystem are also indicated at the top axis. 



DFT. However, for the evaluation of the conductance, 
these studies have neglected the presence of valence d- 
electrons and the scattering due to the non-local pseu- 
dopotential. This approximation is not justified a priori: 
for example, it is clear that the bands due to c?-states 
are very close to the Fermi level in infinite linear chains 
of gold and this indicateSpthat these could play a role, 
especially for a finite bias.EJ 



1. Model 

In this section we consider gold wires connected be- 
tween the (111) planes of two semi- infinite gold elec- 
trodes. In order to keep the computational effort to a 
minimum we will limit our model of the electrode system 
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FIG. 8: Models used for the gold wires calculations. The 
white atoms correspond to contact region C while the grey 
atoms correspond to the L and R regions in Fig. 0. The black 
atoms are only included in the initial Siesta calculation and 
can be added in order to yield a better initial density matrix 
for the subsequent TranSiesta run. We have used 2(A) and 
3(B) surface layers in the contact region. 



to a small unit cell (3 x 3) and use only the F-point in the 
transverse (surface) directions. We have used a singie-C 
plus polarization basis set of 9 orbitals corresponding to 
the 5d and 6(s,p) states of the free atom. In one calcu- 
lation (the wire labelled (c) in Fig. ||) we used double-C 
representation of the 6s state as a check and found no sig- 
nificant change in the results. The range of interaction 
between orbitals is limited by the radii of the atomic or- 
bitals to 5.8 A, corresponding to the 4th nearest neighbor 
in the bulk gold crystal or a range of three consecutive 
layers in the [111] direction. We have checked that the 
bandstructure of bulk gold with this basis set is in good 
agreement with that obtained with more accurate basis 
sets for the occupied and lowest unoccupied bands. 



We have considered two different configurations of our 
calculation cell, shown in Fig. ||. In most calculations 
we include two surface layers in the contact region (C) 
where the electron density matrix is free to relax and 
we have checked that these results do not change signif- 
icantly when three surface layers are included on both 
electrodes. We obtain the initial guess for a density ma- 
trix at zero bias voltage from an initial Siesta calcu- 
lation with normal periodic boundary conditions in the 
transport [z) directionEJ In order to make this density 
matrix as close to the TranSiesta density matrix we 
can include extra layers in the interface between the L 
and R regions (black atoms in Fig. ^) to simulate bulk. 
In the case of two different materials for L and R elec- 
trodes many layers may be needed, but in this case we 
use just one layer. We use the zero bias TranSiesta 
density matrix as a starting point for TranSiesta runs 
with finite bias. 




FIG. 9: We have considered the distances 9.0 (a), 9.3 (b), 
9.6 (c) and 9.9 (d) A between the two (111) surfaces. All 
wires have been relaxed while the surface atoms are kept fixed. 
Distances are shown in A. 



2. Bent wires 

In a r^revious study by Sanchez-Portal and co- 
workers,E3 a zig-zag arrangement of the atoms was found 
to be energetically preferred over a linear structure in the 
case of infinite atomic gold chains, free standing clusters, 
and short wires suspended between two pyramidal tips. 
In general the structure of the wires will be determined 
by the fixed distance between the electrodes and the wires 
will therefore most probably be somewhat compressed or 
stretched. 

Here we have considered wires with a length of three 
atoms and situated between the (111) electrodes with 
different spacing. Initially the wire atoms are relaxed 
at zero voltage bias (until any force is smaller than 
0.02 eV/A) and for fixed electrode atoms. The four re- 
laxed wires for different electrode spacings are shown in 
Fig. 11(a)- (d). The values for bondlength and bondangle 
of the first wire (a), r = 2.57 K, a — 135°, are close to 
the values found in Ref. ^ for the infinite periodic wires 
at the minimum of energy with respect to unit-cell length 
(r = 2.55 A, a = 131°). 

In Fig. |l^ we show the total energy and correspond- 
ing force as evaluated in a standard Siesta calculation 
for the wires as a function of electrode spacing. The 
force just .before the stretched wire breaks has been 
measuredEjO and is found to be 1.5±0.3 nN indepen- 
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FIG. 10; The change in total energy of the relaxed wires 
during the elongation shown in Fig. ^ as calculated in a Siesta 
run. The force determined from the slopes of the line segments 
is shown also. 



dent of chain length. The total transmission resolved in 
energy is shown in Fig. |l^ for zero bias. The conductance 
in units of Go is given by TTot{Ep) which is 0.91, 1.00, 
0.95, and 0.94 for the (a), (b), (c), and (d) structures 
of Fig. ^, respectively. It is striking that the measured 
conductance in general stays quite constant as the wire 
is being stretched. Small dips below 1 Go can be seen, 
which might be due to additional atoms being intrpduced 
into the wire from the electrodes during the puUH It is 
interesting to note that the value for (a) is quite close 
to the conductance dip observed in Ref. |9l| and we spec- 
ulate that this might correspond to the addition of an 
extra atom in the chain which will then attain a zig-zag 
structure which is subsequently stretched out to a linear 
configuration. 

It can be seen from the corresponding eigenchannel 
decomposition in Fig. 12 that the conductance is due 
to a single, highly transmitting channel, in agreement 
with the expflcjjpepts mentioned earlier and previous 
calculations .EJ-EaESEa This channel is composed of the 
L = orbitals.E3 About 0.5 - 1.0 eV below the Fermi 



energy transmission through additional channels is seen. 
These are mainly derived from the = 1 orbitals and 
are degenerate for the wires without a bend, due to the 
rotational symmetry. 

We have done a calculation for a 5 atom long wire. 
The relaxed structure is shown in Fig. |l^. We note that 
while the bondlengths are the same within the wire, there 
is a different bondangle (143° in the middle, 150° at the 
electrodes). We find that the transmission at zero bias 
is even closer to unity compared with the 3 atom case 
and find a conductance of 0.99 Go despite its zig-zag 
structure. 
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FIG. 11; The total transmission of the wires shown in Fig 
vs. electron energy. 



3. Finite bias results 

Most experimental studies of atomic wires have been 
done in the low voltage regime {V < 0.25 V). Impor- 
tant questions about the nonlinear conductance, stabil- 
ity against electromigration and heating effects arises in 
the high voltage regime. It has been found that the sin- 
gle atom gold wires can sustain very large current den- 
sities, wi1|li.an intensity of up to SOjtiA corresponding to 
1 V bias.O Sakai and co-workeraI3il30 have measured 
the conductance distributions (histograms) of commer- 
cial gold relays at room temperature and at 4K and 
found that the prominent 1 Go peak height decreases 
for high biases V > 1.5 V and disappears around 2 V. 
It is also observed that there is no shift in the 1 Go 
peak position which indicates that the nonlinear copduc- 
tance is small. In agreement with this Hansen et alE3 re- 
ported linear current-voltage {I-V) curves in STM-UHV 
experiments and suggested that nonlinearities are related 
with presence of contaminants.. On the theoretical side 
s,p, d-tight-binding calculationsllJO has been performed 
for voltages up to 2.0 V for atomic gold contactSpbetween 
(100), (111) and (110) electrodes. Todorov et aMH has 
addressed the forces and stability of single atom gold 
wires within a single orbital model combined with the 
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FIG. 12: Eigenchannel transmissions of the wires in Fig. ^ 
Only three channels give significant contribution within the 
energy range shown. 
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FIG. 13: Relaxed structure of a 5 atom long chain. Distances 
are shown in A. 
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FIG. 14: The total transmission of all channels and eigenchan- 
nel transmissions of the 5 atom long chain shown in Fig. [l^ . 



fixed atomic charge condition. 

Here we study the influence of such high currents and 
fields on one of the wire structures ((c) in Fig. ^). We 
have performed the calculations for voltages from 0.25 

V to 2.0 V in steps of 0.25 V. In Fig. |l5| we show the 
eigenchannel transmissions for finite applied bias. For 
a bias of 0.5 V we see a behavior similar to the V 
situation except for the disappearance of the resonance 
structure about 0.7 eV above Ep in Fig. 12(c). For 0.5 

V bias the degenerate peak 0.75 eV below Ep which is 
derived from the = 1 orbitals is still intact whereas 
this feature diminishes gradually for higher bias. Thus 
mainly a single channel contributes for finite bias up to 
2 V. It is clear from Fig. |l^ that the transmissions for 
zero volts cannot be used to calculate the conductance 
in the high voltage regime and underlines the need for a 
full self-consistent calculation. 

The calculated I-Vcnvve is shown in Fig. |l6[ We ob- 
serve a significant decrease in the conductance {I /V) for 



high vpltages. This is in agreement with tight-binding 
resultalj where a 30% decrease was observed for a bias,pf 
2 V. For wires attached to (100) and (110) electrodesBO 
a quite linear I-Vwas reported for the same voltage 
range. 

In Fig. |l^ and Fig. |l^ we plot the voltage drop - i.e. 
the change in total potential between the cases of zero 
and finite bias, for the case of 1 V and 2 V, respectively. 
We observe that the potential drop has a tendency to 
concentrate in between the two first atoms in the wire in 
the direction of the current. A qualitatively similar be- 
havior was seen in the-tight-binding results for both (100) 
and (111) electrodeala and it was suggested to be due to 
the details of the electronic structure with a high density 
of states just below the Fermi energy derived from the 
d orbitals (and their hybridization with s orbitals). The 
arguments were based on the atomic charge neutrality as- 
sumption. In the present calculations, this assumption is 
not made, although the selfconsistency and the screening 
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FIG. 15: The eigenchannel transmissions for bias voltages of 
0.5 V, 1 V, 1.5 V, and 2 V for the wire shown in Fig. |. The 
conductance is determined from the average total transmis- 
sion from hl to iir. The voltage window is shown with thick 
dashed lines. 
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FIG. 16: The current-voltage {I-V) curve for the wire shown 
in Fig. |. 
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FIG. 17: The voltage drop for applied bias of 1 V in a plane 
going through the wire atoms. In the surface plot the wire 
atom positions are shown as black spheres. In the contour 
plot below the solid contours (separated by 0.1 eV) show the 
voltage drop. The dashed contours are shown to indicate the 
atomic positions. 

in the metallic wire will drive the electronic distribution 
close to charge neutrality. ThiS|«puld not occur in the 
case of a non-metallic contactED'Ell 

The number of valence electrons on the gold atoms is 
close to 11. There is some excess charge on the wire 
atoms and first surface layers (mainly taken from the 
second surface layers). The behavior of the charge with 
voltage is shown in Fig. The minimum in voltage 
drop around the middle atom for high bias (see Fig. Hq) 
is associated with a decrease in its excess charge for high 
bias. The decrease is found mainly in the s and d^z 
orbitals of the middle atom. 



4. Forces for finite bias 

We end this section by showing the forces acting on 
the three atoms in the wire for finite bias in Fig. We 
evaluate the forces for non-equilibrium inJie same man- 
ner as for equilibrium Siesta calculationsLJ by just using 
the nonequilibrium density matrix and HstHiiltonian ma- 
trix instead of the equilibrium quantities .^3 We find that 
for voltages above 1.5 V that the first bond in the chain 
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FIG. 18: Same as in Fig. |l^ for a bias of 2 V (contours sepa- 
rated by 0.2 eV). 
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FIG. 19: The (Mulliken) change in excess charge (in units 
of the electron charge) on the wire atoms and average excess 
charge on the surface atoms in the first and second left and 
right electrode layers. For high bias the 2nd right electrode 
layer takes up some of the excess charge. 



wants to be elongated while the second bond wants to 
compress. Thus the first bond correspond to a "weak 
spot" as discussed by Todorov et alEXS We note that 
the size of the bias induced forces acting between the 
two first wire atoms at 2 V is-close to the force required 
to break single atom contactstil (1.5±0.3 nN) and the re- 
sult therefore suggests that the contact cannot sustain 



FIG. 20: The forces acting on the wire atoms when the bias is 
applied (the radius of the circles correspond to 0.5 uN). The 
tensile force in the bond between the two first atoms is about 
1 uN for 1 V and 1.5 nN for 2 V. 



a voltage of this_magnitude, in agreement with the re- 
lay experiments.LJ A more detailed calculation including 
the relaxation of the atomic coordinates for finite volt- 
age bias is needed in order to draw more firm conclusions 
about the role played by the nonequilibrium forces on the 
mechanical stability of the atomic gold contacts. We will 
not go further into the analysis of the electronic struc- 
ture and forces for finite bias in the gold wire systems 
at this point, since our aim here is simply to present the 
method and show some of its capabilities. A full report 
of our calculations will be published elsewhere. 



C. Conductance in nanotubes 

Finally, we have applied our approach to the calcula- 
tion of conductance of nanotubes in the presence of poipjb 
defects. In particular the Stone-Wales (SW) defectEj 
(i.e., a pentagon- heptagon double pair) and a monova- 
cancy in a (10, 10) nanotube. The atomic geometries of 
these structures are obtained from a Siesta calculation 
with a 280-atom supercell (7 bulk unit cells), where the 
ionic degrees of freedom are relaxed until any component 
of the forces is smaller than 0.02 eV/A. We use a single-C 
basis set, although some tests were made with a double- 
C basis, producing very similar results. The one dimen- 
sional Brillouin zone is sampled with five k-points. The 
forces do not present any significant variation if the the 
relaxed configurations are embedded into a 440-atom cell, 
where the actual transport calculations are performed. 

In a perfect nanotube two channels, of character tt and 
TT*, contribute each with a quantum of conductance, Gq. 
In Fig. |l| we present our results-fpt zero bias for the SW 
defect. Recent ab initio studiea23'LZl are well reproduced, 
with two well defined reflections induced by defect states. 
The two dips in the conductance correspond to the clo- 
sure of either the n* (below the Fermi level) or the tt 
channel. 
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For the ideal vacancy the two antibonding states asso- 
ciated with broken a bonds He close to the Fermi level. 
The coupling between these states and the tt bands, al- 
though small, suffices to open a small gap in the bulklike 
TT — TT* bands. The vacancy-induced states appear within 
this gap. Otherwise, the three two-coordinated atoms 
have a large penalty in energy and undergo a large re- 
construction towards a split vacancy configuration with 
two pentagons, '-^ 2 eV lower in energy. Two configu- 
rations are possible, depending on the orientation of the 
pentagon pair, depicted in Fig. We have found that 
there is a further 0.4 eV gain in energy by reorienting 
the pentagon-pentagon 60° off the tube axis (Fijg^p2|b) 
resulting in a formation energy of i?/ — 6.75 eV.ES The 
bonding of the tetracoordinated atom is not planar but 
paired with angles of ^ 158°. Some of these structurea 
were discussed in previous tight-binding calculations .l3 
This is at variance with the results of Refs. p9p7| , possi- 
bly due to their use of too small a supcrccU which does 
not accommodate the long range elastic relaxations in- 
duced by these defects. 

The conductance of these defects, calculated at zero 
bias (Fig. does not present any features close to the 
Fermi level. This is in contrast to the ideal vacancy, 
where reflection related to the states mentioned above 
are present. Two dips appear, at possitions similax. to 
those of the SW deffect. An eigenchannel analysisLj of 
the transmission coefficients gives the symmetry of the 
states corresponding to these dips. The metastable con- 
figuration is close to having a mirror plane, containing the 
tube axis, except for the small pairing mode mentioned 
before. The mixing of the tt and tt* bands is rather small. 
The lower and upper dips come from the reflection of the 
almost pure tt* and tt eigenchannels, respectively. This 
behavior is qualitatively similar to the SW defect. On 
the other hand, for the rotated pentagon pair there is no 
mirror plane and the reflected wave does not have a well 
defined character. 



VII. CONCLUSION 

We have described a method and its implementa- 
tion(TRANSiESTA) for calculating the electronic struc- 
ture, electronic transport, and forces acting on the atoms 
at finite voltage bias in atomic scale systems. The 
method deals with the finite voltage in a fully self- 
consistent manner, and treats both the semi-infinite elec- 
trodes and the contact region with the same atomic de- 
tail. 

We have considered carbon wires connected to 
aluminum electrodes where we find good agreement 
with results- published earlier with another method 
(McDCAL)Lfl for electrodes with finite cross section. We 
find that the voltage drop through the wire system de- 
pends on the detailed structure of the electrodes {i.e. 
periodic boundary conditions vs. cross section). 

The conductance of three and five atom long gold wires 
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FIG. 21: Transmission coefficient of pentagon-heptagon dou- 
ble pair vs. energy measured with respect to the Fermi level. 
The dotted line shows the transmission of a perfect nanotube. 

(a) 




FIG. 22: Atomic configurations for the vacancy defect in a 
(10,10) nanotube: (a) metastable and (b) ground state. 



with a bend angle has been calculated. We find that the 
conductance is close to one quantum unit of conductance 
and that this result is quite stable against the bending 
of the wire. These results are in good agreement with 
experimental flndings. For flnite bias we find a non-linear 
conductance in agreement withp^revious semi-empirical 
calculations for (111) electrodes .EJ We find that the forces 
at finite bias are close to the experimental force needed 
to break the gold wiresEll for a bias of 1.5 to 2.0 V. 
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FIG. 23; Similar to Fig. ^ but now for the vacancy: ground 
(continuous line) and metastable (dashed line) configurations. 

Finally we have studied the transport through a 
(10,10) nanotube with a the Stone- Wales defect or with 
a monovacancy (a calculation involving 440 atoms). We 
have found good, agreement with recent ab initio studies 
of these systemsOO 
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